Two-site two-electron Holstein model: 
a perturbation study 



Jayita Chatterjee^ and A. N. Das 



Theoretical Condensed Matter Physics Group 
Saha Institute of Nuclear Physics 
1/AF Bidhannagar, Kolkata 700064, India 



PACS No.71.38. +i, 63.20.kr 



Abstract 



The two-site two-electron Holstein model is studied within a perturbation method 
based on a variational phonon basis obtained through the modified Lang-Firsov 
(MLF) transformation. The ground-state wave function and the energy are found 
out considering up to the seventh and eighth order of perturbation, respectively. The 
convergence of the perturbation corrections to the ground state energy as well as 
to the correlation functions are investigated. The kinetic energy and the correla- 
tion functions involving charge and lattice deformations are studied as a function of 
electron-phonon (e-ph) coupling and electron-electron interactions for different values 
of adiabaticity parameter. The simultaneous effect of the e-ph coupling and Coulomb 
repulsion on the kinetic energy shows interesting features. 



1. Introduction 

The Holstein model [1] is one of the fundamental models describing the interactions 
of conduction electrons with lattice vibrations and has been studied widely for a 
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long time using different analytical and numerical methods. However, conventional 
analytical methods based on variational and perturbation approach are found to be 
much less satisfactory than the numerical methods; particularly in a wide region of 
intermediate values of electron-phonon (e-ph) coupling and hopping. The importance 
of a reliable analytical approach is crucial due to the failure of standard analytical 
approximations on the Holstein model; the Migdal approximation [2] is valid for a 
weakly coupled adiabatic system, while the Lang-Firsov canonical transformation [3] 
could not be suitable far from the nonadiabatic limit and strong coupling limit. The 
failure of the standard analytical techniques in the intermediate range of coupling has 
been observed by many workers [4] while comparing with exact results. Recently we 
[5] developed a pcrturbative expansion using different phonon bases obtained through 
the LF, MLF and MLF with squeezing (MLFS) transformations for a two-site single 
electron Holstein model. We found that the perturbation corrections within the MLF 
and MLFS methods are much smaller than the corresponding LF values for weak 
and intermediate couplings. The MLF method shows a very satisfactory convergence 
in perturbation expansion and is able to produce almost exact results for the entire 
region of the e-ph coupling from weak to intermediate values of hopping. 

Our aim is to extend the applicability of our approach, developed for a single elec- 
tron system, to a many-electron system. To examine that whether this perturbation 
method based on MLF phonon basis works well or not with increasing number of 
electrons and in presence of electron-electron (e-e) interaction, we consider a two-site 
two-electron Holstein model which is the minimal model to investigate the effect of 
e-e interaction on the polaronic properties and the combined effects of the electronic 
correlation and the e-ph coupling on the kinetic energy of the system. It may be men- 
tioned that the polaron crossover has been studied in two-site and four-site systems 
in presence of e-e interaction by different groups [6, 7, 8] with the Holstein model 
previously. But most of the studies were done at the zero phonon approximation 
level with MLF or MLFS transformed Hamiltonian. However, the combined effect 
of e-e interaction and e-ph interaction on the kinetic energy has not been properly 
studied and this is the main motivation of the present study. 

In our study, we have developed a perturbation series using the MLF phonon 
basis and examined the convergence by comparing contributions of different orders of 
perturbation to the energy and the correlation functions for the ground state. The 
variation of different correlation functions involving charge and lattice deformation 
and the kinetic energy of the system with e-ph coupling are presented and discussed 
for different values of adiabaticity parameters and Coulomb repulsion. 
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2. Model Hamiltonian and formalism 

The two-site two-electron Holstein Hamiltonian is 

i,a i,a i 

where i =1 or 2, denotes the sites. Cia (cj^.) is the annihilation (creation) operator for 
the electron with spin a at site i and rii^ (=c|^Cio.) is the corresponding number oper- 
ator, Qi and g2 denote, respectively, the on-site and inter-site e-ph coupling strengths. 
t is the usual hopping integral. 6j is the annihilation (creation) operator for the 
phonons corresponding to interatomic vibrations at site i, ujois the phonon frequency, 
U and V are, respectively, the on-site and inter-site Coulomb repulsion between the 
electrons. 

Introducing new phonon operators a = (&i + b2)/\/2 and d — {bi — b2)/V^ the 
Hamiltonian is separated into two parts (H — Ha + Hg) [5]: 

+ ojog-irii — n2)(d + d}) +u;od^d (2) 

and 

ifa = co'oa'^a - ooon^gl (3) 

where g+ = {gi+g2)/V2, g_ = {gi-g2)/V2, a = a. + ng+, a^ = +ng+, n = ni + n2 
and rii — riii + riii. 

i?a describes a shifted oscillator which couples only with the total number of 
electrons n, which is a constant of motion. represents an effective e-ph system 
where phonons directly couple with the electronic degrees of freedom. We now use 
the MLF transformation where the lattice deformations produced by the electron 
are treated as variational parameters [7, 9, 10]. For the present system the relevant 
transformation is. 

Ha = e^'Hae-'' (4) 

where R = X{ni—n2){d^ — d), A is a variational parameter related to the displacement 
of the d oscillator. 
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The transformed Hamiltonian is then obtained as [7] 
Hd = (^od^'d + J2 ^p'^iT - m AAaC2a exp(2A(d^ - d)) + exp(-2A(d'^ - d))] 

i,a <y 

+ ujQ^g- - A) (ni - 712) ((i + S) + [4 X! ^^T^a + rtxn^ (5) 

i 

where 

ep = e - cc;o(25(_ - A) A (6) 
C/e = f/ - 2cJo(2^- - A) A (7) 
K = y + 2a;o(2y- - A) A (8) 

For the two-site two-electron system there are six electronic states of which three 
belong to the triplet states as follows: 

\i>Tx) = 4T4TlO)e (9) 

|V'T2) = ^ (4^4; + c1^4t) (10) 
IV'T3) = 4i4j0)e (11) 

and other three states are singlet as given below: 

(12) 
(13) 
(14) 

The triplet states do not mix with each other or couple with the d-oscillator. 
In the following we will only consider the singlet states and the rf-oscillator which 
are mutually coupled through the off-diagonal terms of the Hamiltonian (5). For 
a perturbation method it is desirable to use a basis where the major part of the 
Hamiltonian becomes diagonal. For this purpose we choose the basis set |'^+,A/"), 
|'0i,-/V) and iV'-j-ZV) defined as 

|V'+,A^) = ^=i=(a^|V'2)+&iv|V'3))|A^) (15) 
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\^uN) = |V^i)|iV) (16) 
|V'-,A^) = ^=l=(a^|V'3)-&iv|V'2))|Ar) (17) 

where |^/'+), and represent the bonding, nonbonding and antibonding elec- 
tronic states, respectively and |A'") denotes the A^th excited oscillator state in the 
MLF phonon basis where takes integer values from to oo. 



Here a^v = -(f/e - Ve) + ^{Ue - V^f + lQtl{N, N) and hN = 4te(iV, N). 

p\{N-p)\ ' 



where U{N,N) = t^E^L^-l)^^ Nc,, Nc, - ^ 



and te — t exp(— 2A^). 

The hopping term in Eq. (5) has both diagonal and offdiagonal matrix elements in 
the chosen basis (Eqs. 15-17). The states I'i/'e, N) (for e — +,1, — ) are the eigenstates 
of the unperturbed Hamiltonian 

Ho = UJod}d + epUi -^te F{d^d) [c{^C2a + c\^cia\ 

i (y 

[4 X! ^'T^ii + (18) 
% 

It may be noted that t^Fid^d) in the above Eq. (18) is the diagonal part of the oper- 
ator t exp[±2A(d| — dif\ in Eq. (5) in the phonon occupation number representation. 

Since exp \l\[S-d)\ = exp[2Ad^] exp[-2Arf]exp[— 



p\ q 

N 

SO F{dU) = ^(-1)" i^Xf^idUy/ip 
p=i 
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For the choice of the unperturbed Hamiltonian in Eq. (18) the perturbation Hamil- 
tonian {H — Hq) has no diagonal matrix element in our chosen basis. This is not 
achieved when one replaces F{d^d) by 1 and Oat, become A^-independent. The 
convergence in perturbation expansion for the latter case is found to be weaker than 
that when F[d'^d) is included in Hq. 

The unperturbed energies corresponding to the eigenstates |'0e, (e = +, 1, — ) 

are 



E^^N = Nu;o + 2ep + ^^^--^{U,-Ver + ml{N,N), 
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^1,N — 2ep + Ue 

andE% = Nuo + 2e^ + ^ + ^^(£7^ - K)^ + 16t^(iV, iV) (19) 

The part of the Hamiltonian, which is not included in Hq, is taken as the pertur- 
bation Hi. For the chosen basis, the state !'?/'+, 0) has the lowest unperturbed energy 
Ef\. The general off-diagonal matrix elements of Hi between the states \'4'±-, N) and 
l^/^i', M) are calculated as (for (iV - M) > 0) 



{N, ^±\Hi\iJj±, M) = T'iteiN, M) {a'j^b'M + OM^iv) for even (N - M). 

= forodd(N-M). (20) 

(AT, V^ili^ilV^^, M) = -2te{N, M) {d^a^ - b'^b'^) for even (N - M). 

= forodd(N-M). (21) 

{N,i/ji\Hi\ilj+,M) = 2a'j^P{N,M) -2te{N,M) b'^ ior odd (N -M). 

= for even (N-M). (22) 

{N,i/ji\Hi\ilj_,M) = -26'^P(A^,M) -2te(Ar,M) for odd (N-M). 

= for even (N-M). (23) 

1 M 

E[(-i)^ 



where t{N,M) = te(2A)^-^J^ 



(iV-M)! 



^^^''''^ -M(M-l)...(M-i?+l)] 



(A^-M + i?)!i?!' 



P(Ar, M) = cuoig- - X)[VM 5n,m-i + VM+T 5iv,M+i] 
a^v — , = and o^y — 



«Ar + blf) \l {a% + b%) 

Now one has to make a proper choice of A so that the perturbative expansion is 
satisfactorily convergent. In our previous work we obtained good convergence of the 
perturbation series for the value of A which minimizes the unperturbed ground state 
energy. Following the same procedure, i.e. minimizing E^^^q with respect to A we 
obtain 

A = 9-/[l + ..f:, J (24) 
uJoWr-^ + 4:-r\ 

where r = {Ue - K) /2te- 
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The corrected ground state wave function may be written as, 
1 



G 



|V^+,0)+ C],\^Pi.N)+ C^\i^±^N) (25) 

Ar=l,3,.. W=2,4,.. 

where Nq ^ 1+ Y (<^iv)^ + H ( (<^iv)^ + {C^f ) is the normahzation 

iV=l,3,.. Ar=2,4,... 

factor. 

The coefficients C^, and for different values of are determined using 

the matrix elements of Hi connecting different unperturbed states and their energies 
following Ref. [5]. These Coefficients are determined up to the seventh order of 
perturbation while the ground-state energy is found out up to the eighth order of 
perturbation. 

We have also studied the static correlation functions involving charge and lattice 
deformations (niMi)o and (niM2)o, where Ui and U2 represents the lattice deformations 
at sites 1 and 2 respectively. These correlation functions may be written in terms of 
relevant operators as 



(a + a''')^ . , id + S), 
m«i,2 = ni[-2g+ H J ± ni[-A(ni - 712) H J 



The final form of the correlation functions are obtained as 



(26) 



{niUi^2)o = 



2A 



al + hl 



N=l,3,5,.. 



+ 



+ 



E 

i=2,4 

E 



1 



± 



Ng 



+ E 



NCj, 



N=3,5,.. \Ja% + h\ 



{a^Cj^ — b]\fC 
{cin-iCn-i ~ bN-iCj^_i) 



N 



(27) 



The kinetic energy of the system in the ground state Eke = — t < V'gI J2a[c[aC2a 
exp(2A((i^ — d)) -\-c\^Ci„ exp(— 2A((i^ — d))] \iI)g >, has also been calculated using the 
ground state wave function \iPg), evaluated up to the seventh order of perturbation. 
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3. Results and discussions 



Prom the relevant analytical expressions, the quantities of our interest have been 
calculated within the LF and the MLF perturbation methods considering 35 phonon 
states (which is more than sufficient for g- < 2.5 in the transformed phonon basis). 
For all numerical calculations we take e = 0, 512 = (i.e., g+ = g-) in this work. 



In figure 1 we have shown the relative perturbation corrections to the ground-state 
energy, i.e., the ratios of the perturbation corrections of different orders to the un- 
perturbed ground-state energy as a function of g^ for t/uj^ = 1. It appears from 
figure 1 that the convergence is more or less satisfactory in both the MLF and LF 
approaches. Apart from the second order correction, higher order energy corrections 
are very small within the MLF approach except in a small region of g- ~ 0.6 to 0.8 
where the convergence is relatively weaker. The energy corrections of fifth to eighth 
order are, however, negligibly small even in this region. The energy correction of 
any order within the LF method is larger than the corresponding correction within 
the MLF approach. For smaller values olt/ujQ the corrections are much less and the 
convergence is better as expected from our previous study [5]. 

To examine the applicability of our method for higher values of t and also in 
presence of e-e interaction, we plot in figure 2 the relative perturbation corrections of 
different orders to the ground-state energy as a function of g^ for t/uoQ = 2 for both 
zero and nonzero values of e-e interactions. Figure 2(a) shows that the convergence is 
satisfactory and except the second order correction, higher order energy corrections 
are small. In presence of electronic interactions (nonzero U and V) , magnitudes of the 
relative energy corrections increase but the convergence is still good (figure 2(b)). We 
calculate the ground-state energy within the MLF method up to the eighth order of 
perturbation. As pointed out previously, the fifth to eighth order energy corrections 
are so small that the energy calculated up to the eighth order of perturbation may be 
treated as an almost exact result. Wc find that our results exactly match with that 
obtained by the exact diagonalization study [11] for the same set of parameters. 

The on-site and inter-site charge-lattice deformation correlations, (ni'Ui)o and 
{'niU2)o, are determined for the ground state considering up to the seventh order 
of perturbation within the MLF method. In figure 3 we have plotted these quantities 
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Figure 1: Variation of the relative perturbation corrections E^^ /E^^ to the 
ground state energy as a function of the couphng strength (g^) for t/uo = 1.0 
for (a) MLF method and (b) LF method. is the nth order perturbation 
correction to the ground state energy and E^^— E^^q is the unperturbed 
ground state energy. 2,.., 8 denote the value of n in E^\ 
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Figure 2: Variation of the relative perturbation corrections / Eq to the 
ground state energy as a function of the couphng strength {g_) for t/uo = 2.0 
for the MLF method for (a) = and (b) U = 2.0, V = 1.0. E^^'^ is the 
nth order perturbation correction to the ground state energy and E^^— E^^q 
is the unperturbed ground state energy. 2,.., 8 denote the value of n in E^\ 
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against g^. When we examine the convergence by comparing the values of the same 
correlation function, calculated up to different orders of perturbation we find that 
the values obtained considering up to the fourth to the seventh order of perturbation, 
remain almost same. This ensures a good convergence of our perturbation method 
for the correlation function and thus the results obtained are expected to match sat- 
isfactorily with the exact result. The unperturbed LF and the MLF results arc also 
plotted in figure 3 for a comparison, which shows that the unperturbed MLF results 
for the correlation functions are much closer to the exact result than the correspond- 
ing LF result. The relative difference between the exact MLF result and the MLF 
(ZPA) result is much less for {niUi)o than for {niU2)o. This is owing to the fact 
that {niU2)o is a very sensitive function of the perturbation corrections to the wave 
function compared to (niMi)o. The MLF perturbation result predicts a long tail for 
(^11*2)0 in the strongly coupled region where both the LF (ZPA) and the MLF (ZPA) 
give zero value for {niU2)o. This result is a signatTirc of the retardation effect even in 
the strong coupling region where even the MLF (ZPA) do not predict the retardation 
effect. 



To show the crossover from the delocalized (large) to localized (small) polaron we plot 
the quantity X12 (= {{n,i — n2){ui—U2))o/g+) against g- in figure 4(a) where the MLF 
perturbation results for different values of t/cuo are shown. The ZPA results within 
MLF and LF methods are also given for t/wo = 1-0 to compare them with the exact 
perturbation results. The MLF (ZPA) result is very close to the MLF perturbation 
result. The MLF perturbation result shows a smooth crossover from large to small 
polaron as e-ph coupling strength increases. If we identify the crossover point as the 
point of inflection of the curve obtained within the MLF (perturbation) approach we 
find that ior t/uo — 1 the crossover occurs at a critical e-ph coupling g_ (g_ ~ 0.8) 
when U = V = 0. At this point the coefficient of one phonon state (C^^^ in Eq. (25)) 
changes its sign. Similar behaviour has been also observed in case of our previous 
single polaron studies [5]. The critical value of the e-ph coupling where polaron 
crossover occurs increases with increasing value of t/ujo, as expected. 

In figure 4(b) we have shown the effect of e-e interaction on the polaron crossover 
by plotting X12 against g^ for several values of U and V. For this two-site two-electron 
system the polaronic properties we have studied, depend on {U — V) rather than U 
and V independently. From figure 4(b) it is seen that the critical coupling, where the 
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Figure 3: Plot of the correlation functions (a) (niMi)o and (b) {niU2)o with g_ 
for t/cuo — l.O along with the ZPA results within LF and MLF methods. The 
solid curve corresponds to that obtained within the MLF method considering 
up to seventh order correction to the ground-state wave function. 
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Figure 4: The variations of Xi2 = {{ni — n2){ui — U2))o/g+ with g_ within 
MLF perturbation method (a) for different values of t/coo and U — 0; the 
MLF (ZPA) and LF (ZPA) results for t/uo = 1.0 are also shown, (b) for 
different values of U and V ior t/u;o — 1.0. 
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Figure 5: The kinetic energy of the system, obtained using the seventh order 
of perturbation within MLF method, as a function of for different values 
t/ujo (=0.1, 0.5, 1.0 and 1.5). 

polaron crossover takes place, increases and the crossover becomes more abrupt with 
increasing {U — V). Increasing on-site Coulomb repulsion makes the crossover region 
more abrupt and extends the large polaron region. Similar results were also reported 
earher [6] . The inter-site Coulomb repulsion, on the other hand, makes the crossover 
smoother and the crossover takes place at a lower value of e-ph coupling. 



The variation of the kinetic energy of the system with e-ph coupling strength is 
shown in figure 5 for ?7 = V = and for different values of t/ujQ (= 0.1, 0.5, 1.0 
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Figure 6: For t/ujQ = 1.0, the kinetic energy of the system, obtained by the 
MLF perturbation method, as a function of g- for different values of U. 
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and 1.5). We find that our perturbation results for the kinetic energy match exactly 
with the results obtained by exact diagonalization study [11] for the corresponding 
set of parameters. It may be mentioned that for t/ujQ — 1.0, 1.5 the convergence is 
not good in a very small region of around the points indicated by the arrows in 
the figure. Except these narrow regions the convergence of our perturbation method 
is excellent in the rest of the parameter space we studied. 

It is well known that the kinetic energy of a tight-binding electronic system is 
suppressed by the introduction of either e-ph coupling or by Coulomb repulsion. 
An important question one may ask what is the simultaneous role of the Coulomb 
repulsion and the e-ph interaction in lowering the kinetic energy of the system. To 
address this question we have plotted the kinetic energy against e-ph coupling strength 
for different values of U for t/ujo = 1.0 in figure 6. It is seen from figure 6 that for 
large U values there is a flat region where the kinetic energy remains unaffected by 
increase of e-pli coupling. Higher the value of U larger is the flat region. Beyond this 
flat region the kinetic energy is suppressed almost exponentially with increasing g_. 
The occurancc of this flat region and then a rapid fall in the kinetic energy for higher 
values of U may be understood from the results of flgure 4(b) where it is found that 
for larger values of U the polaron crossover is abrupt but there is an extended large 
polaron region (in space). The value of Xi2, which is a measure of the difference in 
the lattice distortions produced at the charge residing site and the neighbouring site, 
is small and almost constant in this large polaron region. As a result a flat behaviour 
in the kinetic energy is obtained. Beyond a critical value of g- an abrupt increase in 
Xi2 would lead to a sharp fall in the kinetic energy, which is evident in figure 6. 



In figure 7 the variation of the kinetic energy with U for different values of g- is 
shown. It is found that for g_ — 1.0 and 1.5 there is a region where the kinetic 

energy increases with U in contrast to the usual behaviour that the kinetic energy 
is suppressed by the Coulomb correlation. The latter behaviour is of course seen 
in major parameter space. To understand this unusual feature of increasing kinetic 
energy with U for finite g-, we plot the variation of Xi2 as well as the kinetic energy 
with U in the inset of figure 7 for = 1.0. It is readily seen that as U increases the 
correlation function Xu decreases which results in an increase in the kinetic energy. 
At the crossover region, Xi2 shows a sharp fall which is refiected by the abrupt rise 
in the kinetic energy. For large values of U, Xi2 is small and the polaron is a large 
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Figure 7: The kinetic energy of the system, obtained by the MLF perturba- 
tion method, as a function of U for different values of for t/uJo = 1.0; The 
inset figure shows the variation of X12/8 with U for g_ = 1.0. 
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Figure 8: The kinetic energy of the system, obtained by the MLF pertur- 
bation method, as a function of U for a particular value of e-ph coupling 
(gf_ = 1.0) for different values of t/cuo (=0.5, 1.0, 1.5 and 2.0). X12/8 is also 
plotted for t/a;o =1.5. 
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one. In this case the lattice distortions uniformly spread throughout the whole lattice 
(hence the charge particle behaves more like a normal electron) and the kinetic energy 
is suppressed with increasing U as commonly expected. 

In figure 8 the variation of the kinetic energy with U for different values of t/cuQ 
is shown for a fixed value of e-ph coupling {g_ = 1.0). The kinetic energy increases 
with U. shows a peak and then decreases. Similar behaviour is also observed in figure 
7. The peak is broader for higher values of t. The initial enhancement of the kinetic 
energy with U is due to rapid spread up of the polaron size or decrease in the value 
of Xi2 whereas in the large polaron region U plays its conventional role of suppressing 
the kinetic energy. With increasing t the kinetic energy increases as well as the large 
polaron region is extended. We have also plotted Xi2 in the same figure for t/ujo = 1.5 
to show the rise in the kinetic energy is related to rapid fall in Xi2- For t/ujo = 2.0, 
the polaron is a large one for g_ — 1 and so the small polaron region is not seen for 
this case. 

In summary we apply the MLF perturbation method on a two-site two-electron 
Holstein model taking the diagonal part of the Hamiltonian (in the momentum space 
of the MLF basis) as the unperturbed Hamiltonian. A good convergence is achieved 
by our approach. The ground state energy and the kinetic energy obtained within 
this method match exactly with those obtained by exact diagonalization study. The 
behaviour of different correlation functions involving charge and the lattice distortions 
are reported as a function of e-ph coupling. Effect of the Coulomb repulsion and the 
adiabaticity parameter on the large to small polaron crossover is studied. On-site 
Coulomb repulsion makes the polaron crossover more abrupt. Our study on the 
kinetic energy in presence of both e-ph interaction and Coulomb repulsion shows that 
for large values of U, the kinetic energy is not suppressed by e-ph coupling in a range 
which is wider larger the value of U. In this range the polaron is a large one with 
almost an uniform value of xu- This makes the kinetic energy insensitive to the e-ph 
coupling. When e-ph coupling crosses a critical value a rapid crossover to a small 
polaron occurs and in this region the e-ph coupling is effective in suppressing the 
kinetic energy. On the other hand for finite e-ph coupling the kinetic energy increases 
initially with U and then decays. This is ascribed to the fact that the former region 
corresponds to a small polaron or a crossover region where the main role of U is 
to spread the size of the polaron and make it more delocalized. This results in an 
increase in the kinetic energy whereas for large U values the polaron is a large one 
and the on-site correlation suppresses the kinetic energy in its usual way. 
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